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ABSTRACT 

The engineering of synthetic gene networks has 
mostly relied on the assembly of few characterized 
regulatory elements using rational design principles. 
It is of outmost importance to analyze the scalability 
and limits of such a design workflow. To analyze the 
design capabilities of libraries of regulatory 
elements, we have developed the first automated 
design approach that combines such elements to 
search the genotype space associated to a given 
phenotypic behavior. Herein, we calculated the 
designability of dynamical functions obtained from 
circuits assembled with a given genetic library. 
By designing circuits working as amplitude filters, 
pulse counters and oscillators, we could infer new 
mechanisms for such behaviors. We also high- 
lighted the hierarchical design and the optimization 
of the interface between devices. We dissected the 
functional diversity of a constrained library and we 
found that even such libraries can provide a rich 
variety of behaviors. We also found that intrinsic 
noise slightly reduces the designability of digital 
circuits, but it increases the designability of oscilla- 
tors. Finally, we analyzed the robust design as 
a strategy to counteract the evolvability and noise 
in gene expression of the engineered circuits within 
a cellular background, obtaining mechanisms for ro- 
bustness through non-linear negative feedback 
loops. 

INTRODUCTION 

Over the past decade, we have witnessed the expansion of 
synthetic biology (1), where the attempts for cell 
reprogramming to perform new tasks have fructified in 



the engineering of several synthetic regulatory circuits 
(2-20). Usually, the design of synthetic circuits has been 
inspired on the use of mathematical models (21,22) and 
empirical engineering rules inferred from natural examples 
(23,24), although requiring in many cases a genetic 
fine-tuning to achieve the desired behavior (25). It is 
expected that the widespread use of libraries of previously 
well-characterized genetic regulatory elements (26-29), 
together with the ability of engineering combinatorially 
those elements (30), will allow avoiding trial-and-error 
procedures, which are not efficient for optimizing and 
implementing complex systems. Those designed circuits 
may be later fine-tuned with directed evolution techniques, 
although there is no general methodology for the de novo 
network engineering. In fact, this bottom-up approach is 
commonly used in other areas of engineering where a set 
of off-the-shelf parts with precise specifications of their 
operating points can be used to engineer sophisticated 
systems, and has been already successful to engineer 
novel biological circuits (12,19). 

Large efforts in generating not only genetic diversity, 
especially libraries of promoters (19,31-36), but also 
post-transcriptional regulatory elements (6,14,37-39) and 
synthetic transcription factors (40,41), encourage to use a 
combinatorial approach to design artificial circuits (see 
Supplementary Data for further details). In addition, the 
quantitative characterization of these regulatory elements 
allows inferring simple phenomenological mathematical 
models, which could be used to construct the model of a 
system that assembles different elements. In that way, 
several synthetic biology-oriented design tools have been 
developed to make available a library of mathematical 
models created from that genetic diversity, together 
with an interface to create gene circuits by wiring 
elements (42-47). Notably, such a genetic diversity is 
translated into a functional diversity when assembling 
circuits, and these circuits could be readily compiled into 
nucleic acid sequences. However, the design is reduced to 
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examine one-by-one all possible combinations (e.g. 
simulating the dynamical behavior), resulting in a 
tedious design process. Thereby, the evolutionary algo- 
rithms and optimization techniques (48-52) allow us to 
automate this process to find the desired circuits and 
finally depict the functional diversity of a library of 
regulatory elements. Our novel approach allows 
assembling models of regulatory elements from a library 
and couples this with an automated design strategy. 

In this work, we tackle for the first time fundamental 
questions that naturally emerge from that approach. What 
functional circuits can we engineer with a given library of 
regulatory elements? What is the diversity of possible 
behaviors and what is the designability (defined as the 
fraction of assembled circuits that follow a given 
behavior) of each one? Is one behavior easier to design 
than others? Certainly, these features depend on the 
employed library. We also wonder what is the sensitivity 
of the results to the regulatory elements; in other words, 
how many functional circuits involve a given regulatory 
element? In addition, we look at the robustness of a circuit 
by locally perturbing its parameters and evaluating the 
resulting fitness. At fixed network topology, we further 
analyze the whole parameter space that provides the 
targeted functionality, which accounts for the robustness 
of all operative points and asymptotically tends to a value 
that we call asymptotic robustness. Indeed, this property 
accounts for the ability to design such a circuit given the 
limitation of the number of genetic elements, and it could 
be important to analyze the natural occurrence of certain 
genetic architectures. All in all, to solve these questions, 
we developed a computational framework to assemble, 
simulate and design circuits, and that allowed us to 
explore the functional diversity that came from the 
assembled circuits with certain behavior (Figure 1). 
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Figure 1. Scheme of the design platform adopted by harnessing a 
library of models of composable regulatory elements. We explore the 
functional networks that can be engineered either by exhaustive 
combinatorial assembly or by heuristic optimization. 



The design of circuits was accomplished by a selection 
step according to a dynamical behavior-based fitness 
function that can also account for robustness. Initially, 
we applied the methodology to design several functional 
circuits with unlimited genetic diversity (given by the par- 
ameter space) and study their asymptotic robustness. 
Then, we designed complex circuits by plugging functional 
modules. Subsequently, we dissected the whole dynamical 
spectrum of a limited library of regulatory elements and 
analyzed the properties of the resulting circuits. We also 
analyzed the dependence of these results on the constitu- 
ent library and how they could change when the 
stochasticity of the cell is taken into account. Later, we 
showed the application of our methodology to design a 
robust circuit. Finally, we discussed the reliability and 
implementability of the designed circuits. 



MATERIALS AND METHODS 

Mathematical model for gene regulatory circuits 

For modeling genetic circuits, we used a coupled system of 
differential equations. We considered three different types 
of species: mRNA (it can also be non-coding), proteins 
(mainly transcription factors) and small molecules that 
interact with proteins to activate or inhibit their regulatory 
ability (e.g. isopropyl (3-D-l-thiogalactopyranoside — 
IPTG — inhibits the activity of LacI). Likewise, the 
production of the i th mRNA from a regulated 
promoter follows 



dxj 
lit 



CjiypUj) - (8+n)x h 



(1) 



where the term f(yj,Uj) is the transcription rate (yj and uj 
represent the concentrations of the /-th protein and its 
regulating chemical, respectively), C the gene copy 
number, S the mRNA degradation coefficient and \i the 
growth rate of the cell (dilution term). C = 1 is assumed to 
be constant in this work. For the computational design of 
a circuit, we did not impose variability on \i but 
we assumed a constant value (e.g. ii = 0.02 mm -1 ). For 
simplicity, we assumed that all genes in an operon (i.e. 
controlled by the same promoter) have the same mRNA 
expression. The term f(ypUj) accounts for protein-DNA 
and protein-molecule interactions (22), and for constitu- 
tive promoters it is constant (see Supplementary Data for 
further details). Importantly, our approach is independent 
of the choice of this function, thus giving a big degree of 
freedom to the kinetic characterization from experimental 
data. Afterwards, the production of /-th protein (y,) is 
given by 



dt 



= g(Xi,Xj) - (yS+/z)ji 



(2) 



where the term g(x h Xj) is the translation rate and p the 
protein degradation coefficient. The term g(x i: Xj) 
accounts for post-transcriptional regulatory mechanisms, 
such as riboswitches, allowing a further genetic element, 
such as a (ranj-RNA, to control translation (6). In case 
of no post-transcriptional elements, the translation 
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rate is proportional to the mRNA concentration. In 
addition, in this work we only considered first-order 
degradation kinetics. The description of the construction 
of the stochastic model is detailed in the Supplementary 
Data. 

Library of models of regulatory elements 

To construct a library, each genetic regulatory element 
was modeled by transfer functions that related the 
output to the input values. These functions can be fitted 
from experimental data. As DNA fragments, mathemat- 
ical models can be assembled in a standard way to 
simulate the behavior of circuits. Here, we only allowed 
joining promoter and genes, or genes and genes (i.e. two 
consecutive promoters was not allowed; such a construc- 
tion should be specified as a whole part). One useful 
format to store a mathematical model (molecular 
species, kinetic parameters and DNA sequence) is 
systems biology markup language (SBML) (53). Hence, 
we had a single SBML file for the model of each biological 
part; similarly as crystallographic data is stored in protein 
data bank (PDB) format. The models for promoter parts 
only account for the transcription rate, whereas for gene 
parts the model accounts for the translation rate and the 
degradation and dilution rates of mRNA and proteins. 
We selected Hill-function models because their over- 
whelming use in current characterization of transcription 
regulation works. In the future, when more advanced 
models may be used to fit characterization data, 
they could be readily used with our computational 
design procedure. A range of variation can be specified 
for some kinetic parameters; likewise, the corresponding 
value will be susceptible to be changed during the design 
process. 

Exhaustive versus heuristic design 

Multiple circuits can be constructed by harnessing the 
available regulatory elements. To computationally 
explore the functional diversity that offers such a library 
and the designability of certain behaviors, two different 
strategies can be adopted, and our approach provides an 
automated implementation of them. On the one hand, 
following the exhaustive design strategy, all possible 
circuits, up to a maximal number of elements, are 
constructed and simulated. Having the large collection 
of dynamics, a post-processing step is applied to find 
those circuits that behave according to the design specifi- 
cations. This approach allows obtaining the whole 
functional diversity and designability. On the other 
hand, a heuristic design strategy provides a probabilistic 
sampling frame of the functional diversity. It allows itera- 
tively assembling models of existing elements and 
evaluating the performance of the resulting circuit accord- 
ing to a dynamical behavior-based fitness function (54). 
For that, we used Monte Carlo Simulated Annealing 
(MCSA) as optimization scheme (55). A movement in 
the fitness landscape consists in a replacement, addition 
or deletion of a given regulatory element. To evaluate the 
fitness function, we first calculated the average distance 
(metric function) for all genes i considered as outputs, 



for a given target behavior k, between the current circuit 
dynamics (y ik ) and the target one (z lVt ) according to 



fo |logQfr(0) - log(za(0)|xtt(Q<ft 
foX,k(t)dt 



(3) 



where T is the final time (e.g. the time to reach the steady 
state). The metric is in logarithmic scale to properly 
balance the species concentrations, since they can vary 
in several orders of magnitude in biological systems. The 
function x(t) is a weighting factor to only evaluate the 
circuit dynamics in a specified temporal domain (x(t): 
[0,7] — > [0,1]). Subsequently, the fitness function that 
aggregates all targets we used reads 



00 



(4) 



where </> 0 is a normalization constant to adjust the fitness 
value to the metric function (e.g. <p 0 = 3), and Yk gives the 
scalar weight in logarithmic scale for optimizing target k 
(e.g. Yk = 10y/ indicates that target k has 10 times more 
priority than I). If 4>k>4>o f° r one k then we assumed 
xjr = 0. Importantly, this fitness function (v|/ belongs to 
the interval [0,1]) penalizes those circuits that do not 
satisfy simultaneously all targets. Being Ajjr the fitness 
update after a movement, this is accepted with probability 
max{l, exp(z1i/^/7 , McsA)} 5 where T MC sa is the MCSA 
temperature. 7~mcsa is continuously adjusted during the 
optimization process following an exponential cooling 
scheme. 



RESULTS AND DISCUSSION 

Circuit design and asymptotic robustness 

Initially, we constructed a library of artificial regulatory 
elements, including all types of logic combinatorial pro- 
moters of two entries. Additionally, the kinetic parameters 
characterizing those elements were specified as a range of 
variation. This feature allows that the genetic sequence of 
many biological parts could be easily modified to create 
a new part with diminished binding affinity or stability by 
a single mutation. Otherwise, it is much more difficult to 
find a suitable mutation that would increase the binding 
or stability. Therefore, by allowing this range in the 
parameter space, we would enlarge the search space 
while still maintaining the linking with the genotype, 
because the parts from an optimal solution could be 
readily engineered to follow a model agreeing with the 
designed parameters. The nominal values were taken 
from several experimental studies (11,12,27,36). Thus, 
the genetic diversity was almost unlimited, being the 
design space defined by topological and parameter modi- 
fications of the circuit. To explore this space we applied 
the heuristic design strategy to find the optimal assemblies 
of elements and parameterizations that gave functional 
circuits. First, we repeatedly applied the optimization 
method to design all possible circuits relying on a 
feed-forward loop (FFL) structure for one-stripe pattern 
formation (56—59). We found six different architectures 
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Figure 2. Schemes of several FFL-based gene circuits optimized to operate as amplitude filters (also called band detectors). The mathematical model 
for each circuit is provided in SBML format in the Supplementary File sbml.zip. 



for working as an amplitude filter (Figure 2), where five of 
them corresponded to incoherent FFLs (all except 2A), 
and in architectures 2B and 2E repression dominated 
over activation. Moreover, in Supplementary Figure SI, 
we show the circuits corresponding to inverse amplitude 
filters. Certainly, the two regulatory branches with 
opposite sign are responsible for such a behavior, and 
the combinatorial promoter of the downstream 
promoter is central to get a variety of functionally analo- 
gous circuits. Interestingly, some of those architectures 
have been found involved in developmental processes 
(60,61). In Supplementary Figure S2, we illustrate a 
possible implementation of a FFL-based circuit respond- 
ing to intermediate concentrations of IPTG, together 
with its characteristic transfer function. Notably, we 
did not exhaustively construct all possible FFL circuits 
from the library for their scoring. Instead, we prob- 
abilistically sampled the fitness landscape and we 
always found a solution corresponding to one of the 
six FFL structures presented. Moreover, this approach 
can be applied to design functional circuits without 
accounting for the designability of the desired behavior, 
and then study the intrinsic properties of the circuit 
irrespective to the library, such as its asymptotic 
robustness. 

We then investigated the asymptotic robustness of those 
FFL circuits, which functioned as amplitude filters with 
a fold-change (F) of at least one order of magnitude at the 
detection point (F> 10). By constraining the sign of the 
regulations (fixed topology), we obtained a parameter 
space of ~2xl0 5 different combinations for each 
topology shown in Figure 2. Accordingly, the highest 
asymptotic robustness was reached by the architecture 
2E with the 21.29%, followed by the architectures 2A 
with the 17.17% and 2F with the 17.66%, indicating 
that those circuits have a one-stripe pattern-prone 
structure. Interestingly, this could be because the input 
gene (X) has a non-monochromatic regulatory mode 
(i.e. both activator and repressor) in these topologies. 



On the contrary, the architecture 2C was highly sensitive 
to parameter variations with asymptotic robustness of 
only the 2.54%. The architectures 2B and 2D with the 
8.60% and 7.69%, respectively, were in between. 
However, despite of its low asymptotic robustness, the 
structural core 2C is broadly found in many natural 
systems. For instance, in the Drosophila patterning circuit- 
ry, gene hb represses both genes kni and Kr and kni 
also represses Kr (60). In addition, the core 2B is the 
FFL motif most abundant within the regulatory map of 
bacteria and yeast (62). That the genes involved in the 
structures 2B and 2C have a monochromatic regulatory 
mode could explain the increasing presence of these 
circuits. Moreover, from a synthetic perspective, pro- 
moters type NOR and IMPLIES could be engineered by 
placing contiguously the corresponding operators in the 
promoter region (7,12,36). 

Subsequently, we used the optimization method to 
design a circuit able to count. This is an interesting 
example that could already unveil many of the issues we 
meet in more complex networks. Cells may take advantage 
of this sort of circuits to regulate fundamental processes, 
such us telomere length control (63), where a machinery to 
count molecules or events is required. Counters have 
different stable states and rely on memory-like architec- 
tures that allow retaining the initial state, unless a perturb- 
ation switches the system (3,14,18,64). Generally, the 
underlying mechanism of biological counters consists 
in overcome certain threshold after a specific number of 
consecutive pulse-like events. Herein, we attempted the 
design of a two-pulse counter, where we imposed that 
the system had to reach three different states. We 
applied the optimization method to design all possible 
two-gene circuits. We found that, within a delimited 
time domain, all possible circuits were functional and 
reached three states. However, those circuits based on 
an activator-repressor core had a meta-stable state, 
which falls into the basin of attraction of one of the two 
stable states after certain time (Supplementary Figure S3). 
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In Figure 3, we show the phase diagrams for the circuits 
based on a monochromatic regulatory core and one 
self-activation (see in Supplementary Figure S4 the 
corresponding phase diagrams for the circuits based on 
an activator-repressor core). The addition of another 
self-activation on the buffer gene allows having a symmet- 
ric multi-stable device (64). We then computed the asymp- 
totic robustness of these two circuits, by exploring 
exhaustively all possible parameterizations (for that we 
discretized the parameter space into ~2 x 10 5 different 
combinations). The double repression core allowed 
tristability in the 0.2422% of the cases, whereas the 
double activation core in the 0.1787% (relatively low in 
both cases). 

Next, we attempted the design of a two-pulse counter 
relying on just two states. As design specifications, we 
imposed pulses of lOmin within an interval of 50min 
with amplitude of 100-fold. The designed circuit with a 
plausible implementation counting pulses of IPTG is 
shown in Supplementary Figure S5. However, the func- 
tioning of this system (i.e. number of pulses it is able 
to count) depends on the pulse length and interval. 
In addition, we attempted the automated design of a 
tunable genetic timer. These devices consist of memories 
that change the state of operation according to an external 
signal and the time to accomplish this transition (time to 
reach the steady state) can be modulated by another signal 
(19). The designed circuit with a plausible implementation 
is shown in Supplementary Figure S6, together with its 
characteristic transfer function. This circuit consisted in 
a coherent FFL coupled to a memory-like mechanism 
based on a self-activation, and it existed a threshold 
in the IPTG concentration from which the circuit 
responded to different levels of it. 

Hierarchical design and modularity 

Once a functional genetic device is obtained, either from 
computational or rational design methods, and experi- 
mentally validated, it can be integrated in the library as 
a new element to be used in the construction of more 
complex systems. As a first approach, we included in 
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Figure 3. Schemes of two two-gene circuits designed to reach 
tristability, showing the corresponding phase diagrams. Filled and 
open circles represent stable and unstable states respectively. The math- 
ematical model for each circuit is provided in SBML format in the 
Supplementary File sbml.zip. 



our library of regulatory elements a circuit previously 
optimized to operate as a tristable. Remarkably, the 
incorporation into the design procedure of black-box 
modules enhances the optimization of the impedance 
matching, where the output of a device serves directly as 
the input of a downstream one, and could considerably 
enlarge the functional diversity of the library. Hence, 
following such a design approach, we were able to 
obtain complex functions with modular systems. In our 
particular case, we designed a system coupling a tristable, 
an amplitude filter, and a frequency-tunable oscillator 
(Figure 4). Initially, this tristable gave a low concentra- 
tion. After a pulse of 20min with amplitude of 1000-fold 
in the inducer, the device switched its state to reach the 
intermediate concentration level, and subsequently the 
oscillator changed its frequency and the amplitude filter, 
which operates as a detector of the intermediate state, 
reached its ON state. After a second pulse, the device 
switched to its high concentration point, inducing a new 
change in the frequency of the oscillator and giving the 
detector back to its OFF state. Interestingly, the 
frequency-tunable oscillator evolved to couple two differ- 
ent regulatory mechanisms, and the external signal 
switched from one mechanism to another, then changing 
the frequency of the oscillations. In addition, with the 
consideration of delayed reactions (e.g. due to translation 
and multimerization) we could obtain complex oscilla- 
tions, which can drive to a route toward chaos (65). 
In fact, this mechanism has been previously applied to 
design genetic oscillators with a minimal number of 
elements (16). In Supplementary Figure S7, we show 
a genetic clock designed by optimization coupled to a 
tristable, which modulates its frequency and the shape of 
its oscillations according to the states of the switch-like 
circuit. 

However, one important issue in such an approach is 
the possibility of the loss of function of a device when 
plugging it to a downstream module. This effect, usually 
called retroactivity (66), emerges when a transcription 
factor plays two different roles in both modules, and 
is indeed a consequence of the limited protein amount. 
This result may have significant consequences on the 
dynamics of the system, even when the stochasticity of 
the cell is taken into account (67). Here, our modeling 
neglects this effect by assuming that the concentration of 
free protein is always much higher than the protein bound 
to DNA (22); also as an imposition to ensure modularity 
in the design and to be able to combine different elements 
from the library. Although for many systems this 
approach is valid (9), it could be found some examples 
where such a model is not too accurate. To solve this 
problem in practice, one strategy would be to impose as 
a design constraint that the output gene had no regulatory 
effects on the circuit. Likewise, this output could be used 
as the input in further downstream modules with increased 
guarantees of a proper functioning. Thereby, in our 
system of Figure 4, gene U could be split into two genes, 
one for working within the tristable device and another 
for setting the amplitude filter and the oscillator, although 
still it would exist a coupling between these two devices 
due to a common input. 
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Figure 4. (A) Scheme of a complex regulatory system comprising a frequency-tunable oscillator and a state detector, designed by using the tristable 
device as an element of the library. Moreover, we show the transfer functions of the different devices that form the system. (B) Dynamics of the 
output genes of the complex system. Pulses in the input (/) of 20min and 1000-fold of amplitude were applied at t - 1000 and 2000 min. 
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Functional diversity and designability 

We further studied the designability of a given dynamical 
behavior. For that, we constructed a library of SBML 
models of well-characterized regulatory elements previ- 
ously implemented in vivo. Likewise, the corresponding 
kinetic parameters were fitted from experimental data 
and kept fixed. Using this library, we constructed by 
in silico assembly all possible architectures up to three 
genes, giving 501 952 different circuits (see the different 
configurations in Figure 5A). We systematically imposed 
X-cl as output gene in all circuits. Thereby, we computed 
the dynamics of all circuits to perform an analysis of the 
behaviors that could be obtained with such a library. For 
this work, we considered a library of 36 elements, 
involving 5 genes and 31 synthetic promoters. As genes, 
we contemplated the classical repressors LacI, TetR and 
A.-cI, and the activators AraC and LuxR. Moreover, we 
built a library with three constitutive promoters with dif- 
ferent transcription rates, 16 single promoters involving 4 
lacO, 4 tetO, 2 araO, 2 luxO (36), 2 X r O (1 1) and 2 A. rm O 
(5) and 12 combinatorial promoters involving 4 lacO- 
tetO, 3 araO-lacO, 2 araO-tetO (36), 1 A RM 0-lacO (12), 
1 1uxO-A r O (7) and 1 luxO-lacO (68). The models also 
accounted for the external molecules (IPTG, anhydrote- 
tracycline -aTc- L(+)-arabinose and acyl homoserine 
lactone — AHL — ) that modify the regulatory ability of 
the transcription factors and represent the inputs of 
the circuits. These SBML models are provided as 
Supplementary Data and they contain the corresponding 
kinetic parameters. Then, for each external inducer we 
considered three different states (low, intermediate and 
high), giving 81 environmental conditions for all combin- 
ations, and four more conditions in which the inducers 
had a pulse-like dynamics. 

By compiling all numerical results (details in 
Supplementary Data), we were able to dissect the dynam- 
ical spectrum of the library (i.e. its functional diversity), 
which included circuits operating as oscillators, amplitude 
filters, memories and different logic gates (Figure 5B). 
As expected, the majority of the circuits functioned as 
logic gates, and because the external signals (IPTG, aTc, 
arabinose and AHL) always activated transcription, the 
set of NAND and NOR gates was highly reduced. In 
addition, ~1% of the assembled circuits was able to 
exhibit oscillations. Furthermore, we found amplitude 
filters in the 0.016% of the cases, and memories in the 
0.436%. Certainly, this spectrum depends on the value 
of F specified to differentiate between two concentration 
levels (F gives their ratio). Herein, we imposed F> 10, 
although we also performed a screening to see the effect 
of different values of F. Not surprisingly, as higher is 
F, the number of functional circuits decreases 
(Supplementary Table SI). In addition, we studied the 
effect of the initial condition on the output gene finding 
that the results were almost independent of this. Certainly, 
the initial condition only affects in memory-like circuits, 
but this effect was captured by imposing pulse-like 
dynamics on the input. Interestingly, the repertoire of 
designed circuits was essentially based on minimal cores 
that provided the required functional mechanism 




Figure 5. (A) Graphical representation of the exhaustive design 
strategy. Starting from a library of composable genetic regulatory 
elements (mathematical models provided in SBML format in the 
Supplementary File sbml.zip), we constructed all possible circuits up 
to three genes for simulation. (B) Dynamical spectrum of the library by 
exhaustive exploration (functional diversity). We represent the percent- 
age of circuits that behave as oscillators, amplitude filters, memories 
and logic gates (designability). To differentiate between two states of a 
circuit, we imposed at least one order of magnitude in concentration. 



(Supplementary Figure S8). These cores illustrate the 
design principles in which the dynamical spectrum is 
based on. However, the use of a limited library and a 
partial set of input conditions, while allowing an exhaust- 
ive exploration, prevent obtaining a comprehensive 
analysis of the design principles. For instance, as we 
have shown above, a double activation core gives a 
memory-like mechanism but it was not found in the 
repertoire of circuits. In addition, all amplitude filters 
were based on the FFL architecture 2C, although 
further circuits, not necessarily FFLs, can be employed 
to read morphogen gradients (59). We did not obtain 
further topologies because the monochromatic regulatory 
mode and the lack of cooperation between transcription 
factors. 

Scalability and sensitivity of designability 

Furthermore, we investigated the dependence of the 
designability of a function on the existing elements of 
the library by calculating the degree of sensitivity of 
each regulator over the resulting dynamical spectrum 
(Figure 6). Accordingly, LacI appeared to be the most 
important regulator, indeed for this particular case of 
study, since it participated in the majority of the 
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Figure 6. Sensitivity analysis of the dynamical spectrum. We release 
one regulatory element of the library (in particular, one gene) to 
analyze its contribution to the dynamical spectrum (we represent the 
remaining number of functional circuits relative to the total). 



functional circuits. In the specific case of the amplitude 
filters, since their mechanism relied on two different re- 
pressions (Supplementary Figure S8), Lad and TetR 
participated in all circuits. Certainly, the addition of 
more regulatory elements in the library would enlarge 
the designability of the different behaviors, and the iden- 
tification of the regulatory cores in Supplementary Figure 
S8 would lead to rationally decide on the elements of more 
interest. In addition, we studied whether the designability 
could be estimated by sampling a small subset of 
assembled circuits instead of an exhaustive exploration. 
This would provide further support to the heuristic explor- 
ation by means of optimization methods. Interestingly, we 
found similar results for the dynamical spectrum of the 
library when analyzing the dynamics of about 1000 
circuits (corresponding to the 0.2% of the total circuits) 
(Supplementary Figure S9). This suggests that even a 
small fraction of assembled circuits is representative of 
the whole population of circuits. By exploiting this fact, 
we could analyze the functional diversity and designability 
of several libraries of models at a minimal computational 
cost or we could study how to enrich the library with 
new regulatory elements. 

We further studied the designability of the different 
behaviors when considering the stochasticity inherent to 
the cellular processes (we focused on intrinsic noise, details 
in Supplementary Data) (69). Since the stochastic simula- 
tion entails a higher computational cost, we considered 
a subset of circuits as before to perform this study, and 
because it is expected this will not strongly affect the 
results. For each condition of inputs, we considered 
the average value and standard deviation of the output 
(computed using the time dynamics after a transient 
period). In general, we found similar results as in the 
deterministic regime (Supplementary Figure S10). We 
could explain this by the fact that in most cases gene 
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expression is high enough, which minimizes the effect of 
intrinsic noise, although in some cases a particular circuit 
topology could also help in such a noise reduction (59). 
However, we found an increase of almost a doubling in 
the number of oscillators. By examining the circuits, we 
realized that circuits based on an activation-repression 
mechanism with fast damped oscillations in the determin- 
istic regime and that were identified as stable circuits were 
then selected as noise-induced oscillators. For the other be- 
haviors, the designability results in the stochastic regime 
were slightly lower. The maximal reduction in 
designability was of ~20% in the case of YES/NOT 
gates. In the circuits that were selected according to the 
deterministic solution but not to the stochastic one, there 
is an increase of noise in protein expression that prevents 
identifying different states of operation. 

Afterwards, we wondered whether a unique circuit 
could exhibit different behaviors. Interestingly, we 
observed special circuits that displayed multifunctionality 
according to different input conditions (e.g. oscillators 
working as amplitude filters, memories or logic gates). 
Supplementary Table S2 shows the number of circuits 
with two functions and the corresponding statistical 
significance. For instance, the 0.3% of the total set of 
circuits functioning as oscillators and memories held the 
two functions by properly setting the environmental 
factors (statistical significance assessed by bootstrapping). 
In addition, we calculated the number of circuits with 
multifunctionality (Supplementary Table S3) showing a 
tendency log-normal in the distribution (slope = —1.463, 
R 2 = 0.982, P<0.1). This sort of circuits is appealing for 
cellular regulation and organization because the circuitry 
rewiring instrumental to change the function is accom- 
plished by an on-the-fly reprogramming sentence (70). 
As well as a single gene can attain several functions 
[e.g. a protein with different enzymatic properties (71)], 
a multifunctional circuit can be exploited by the cell to 
exert a conditional control of different responses. 

Robust design for a cell environment 

One of the main problems in synthetic biology is the 
reliability of the engineered circuits when they are 
deployed in a continuous cellular background. Certainly, 
if there is not a selective pressure to maintain a circuit and 
it results in a cost for the cell, the genetic parts that form it 
will accumulate mutations (involving point mutations, in- 
sertions and deletions) and a quantity of them will entail 
the loss of function of the circuit (27,72). Thereby, 
we implemented a strategy based on robust design to 
counteract this undesirable effect by which the resulting 
circuit has the ability of maintaining the same behavior 
under perturbations in the kinetic parameters of the model 
(i.e. a sort of mutational robustness) (54)]. This also deals 
with the eventual parameter uncertainty of the inferred 
models. To this end, we expanded the metric function to 
account for robustness, scalarizing the corresponding 
multi-objective optimization problem, given by 

(Pt = (l-X)ct> k +X(<t> k ), (5) 
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where X is the degree of robustness specified (i.e. for X = 0 
no robustness is imposed), and (4>k) is the average of the 
metric function over a large number of perturbations. As a 
case of study, and since the promoter region is usually a 
preferential site to accumulate mutations (72), we just 
focused on perturbations in the binding affinities that 
randomly change the parameter values up to 2-fold 
(changes in all parameters simultaneously). In addition, 
to accelerate the convergence future work could just 
evaluate that function in a more reduced neighborhood 
with a first-order approximation (73). 

Accordingly, we attempted the design of a robust amp- 
litude filter satisfying the (input, output) specifications of 
{(1,0) a (100,100) a (10 4 ,0)} (Figure 7A). The architecture 
of the resulting circuit involves a FFL and for its imple- 
mentation we would need a promoter type NAND, which 
could be engineered by using artificial operator sites that 
would be recognized by heterodimers from chimeric 
proteins (74). We then used this circuit as starting point 
(keeping fixed the topology) to further apply the optimiza- 
tion procedure with X = 0. As we could observe, the new 
circuit reaches a higher fitness at X = 0, which indicates 
that it exists a cost for robustness (Figure 7B). While being 




a, (robustness weight) 

Figure 7. (A) Scheme of a genetic circuit optimized to operate as a 
robust amplitude filter. U and Z are the input and output, respectively. 
(B) Fitness function versus the robustness weight for the robust circuit 
shown above (red) and the same circuit optimized with A = 0 (blue), 
illustrating the cost of robustness. 



this cost low, this allows the circuit to maintain an 
operative mode even under genetic perturbations. 

In addition to the evolvability of living systems, gene 
circuits are subjected to stochastic processes that may 
prevent the expected behavior. The use of stochastic 
models could allow designing more robust circuits, or at 
least circuits that maintain their deterministic behavior. 
One way to explore the effect of noise in the design of 
circuits would be to in silico evolve a circuit to be robust 
against noise and compare it with a circuit not evolved 
with such a selective pressure. We applied our optimiza- 
tion method to design circuits that maintained a targeted 
dynamics in presence of noise (both intrinsic and extrin- 
sic), obtaining different circuits with self-repressions and 
low translation rates (Supplementary Figure Sll) (69). 
On the contrary, this approach could be applied to 
design circuits displaying high noise levels (75), since 
they may provide to the cell certain advantage under 
unpredictable environmental conditions (20). 



CONCLUSION 

In this work, we have tackled the problem of the 
designability of a given gene dynamics provided a 
library of composable regulatory elements, considering 
that the functional circuits come from combining different 
elements of the library. This measure of designability 
quantifies the entropy of a given dynamical behavior. 
For that, we have developed a computational method- 
ology that allows exploring the diversity of behaviors 
that can be obtained by assembling circuits by means of 
two different design strategies: one based on heuristic op- 
timization and other based on exhaustive simulation of 
circuits. We have taken advantage of current characteriza- 
tions of regulatory elements into libraries of mathematical 
models (26-29), allowing to rapidly select the regulatory 
element of interest for our circuit. Although the emergence 
of unexpected behaviors is always an issue in synthetic 
biology, it is anticipated that the use of standardized 
parts allows reducing the endless tweaking process when 
engineering a synthetic gene circuit (12,19). Using a proper 
mathematical formulation, we were able to generate a 
large collection of genetic circuits by assembling those 
regulatory elements, and identify the functional subset ac- 
cording to certain specifications. Initially, we constructed 
an artificial library of models to design circuits by opti- 
mization toward a configuration satisfying the specifica- 
tions. We designed filters and counters of gene expression, 
which allowed us to find new regulatory mechanisms 
able to provide such behaviors. Sometimes the behavior 
requires a very precise genotype, making improvable to 
get many cells with such behavior in a heterogeneous 
population. To investigate this, we have defined the 
concept of asymptotic robustness, which provides a 
measure of the maximum genotypic heterogeneity for 
a given of phenotypic behavior. In the long term, it is 
expected that synthetic biology projects will provide 
many examples of standardized circuits with a given 
dynamics, which could be incorporated into the available 
libraries. Then, one could extend our analysis to such 
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cases. One issue here would involve the interfacing of such 
modules. To analyze this, we exploited one of our 
designed circuits as a single element of the library to 
obtain a complex system involving such a functional 
unit, illustrating a hierarchical design approach and 
allowing the design of plug-and-play devices with 
optimal impedance matching. 

Given a library of regulatory elements, it is possible to 
construct many circuits with various dynamical behaviors. 
But some behaviors occur more often than others. To 
quantitatively analyze this, we computed the designability 
of a set of useful behaviors. There, we constructed a more 
reduced library of regulatory elements to assemble all 
possible circuits up to three genes and process their 
dynamics. Remarkably, the library involved promoters 
that had been previously characterized and even used for 
engineering various synthetic circuits in the bacterium 
Escherichia coli. Interestingly, we found that a limited 
library could encode a large number of behaviors. 
Certainly, our computational method allowed construct- 
ing and simulating the dynamics of this large set of circuits 
and assisted to dissect the spectrum of dynamical behav- 
iors and study their designability. We found that a same 
genotype could have several functions depending on the 
external signals. Nevertheless, as the size of the circuits 
and the number of elements of the library increases, the 
exhaustive design strategy becomes unpractical, thus 
requiring heuristic methods. Since noise is an important 
factor that affects the dynamics of a circuit, we also 
included it in our analysis. The consideration of intrinsic 
noise slightly reduces the designability of digital circuits, 
but it increases the designability of oscillators. This is 
understandable from the fact that digital devices are 
steady-state based circuits, where noise could only spoil 
the behavior. On the other hand, oscillatory circuits are 
dynamical systems, where the noise could contribute to 
enhance the behavior. In addition, we expect that our 
results would be maintained when the library is 
enhanced by incorporating more accurate experimental 
measurements of the transcription regulation elements. 
That the new models could be more elaborated and 
accurate, they would not much change the fitness land- 
scape and thus the designability of behaviors. 

Herein, as opposite to other optimization-based design 
approaches (48-52) where it is difficult or even impossible 
to compile into a reliable DNA sequence the designed 
circuit, our automated approach scrutinizes in a combina- 
torial way the available genetic diversity to optimize 
a circuit with the desired function, and finally output 
a DNA sequence that encodes this designed circuit. 
However, the evolvability of genetic circuits in a continu- 
ous cellular background is a delicate issue. Importantly, 
the robust design approach entails that the expected 
behavior can be also observed throughout the cells of a 
population. Indeed, natural systems are robust to coun- 
teract mutations that, in most cases, entail unpredictable 
changes in the function of these elements (76). Moreover, 
as higher is the load of the circuit for the cell (i.e. quantity 
of resources required for function expression), faster will 
be the loss of function of the circuit because deleterious 
mutations will be selected (27,72). Hence, we tackled the 



design of robust circuits by accounting for parameter 
sensitivity into the objective function. We illustrated that 
such a circuit pays a cost by decreasing its fitness in order 
to gain robustness, which will lead to a higher evolution- 
ary stability. Alternatively, or even in combination, we 
could design circuits that could, while preserving its 
central function, induce a benefit to the cell, by accounting 
for the interface between the circuit and the cellular 
interactome (77), and then assess the selection and 
reliability of the engineered circuit. 

Interestingly, one possible extension to our work would 
be the development of a more complex, hierarchically 
distributed design platform (26). Herein, more diverse, 
characterized regulatory elements would be considered, 
involving transcriptional, riboregulatory, metabolic and 
signaling elements. These different regulatory elements 
would be combined to yield complex functional genetic 
circuits, involving different regulatory mechanisms. 
In addition to new elements, inherent effects such as the 
variation of cell growth rate due to different culture 
media, the delay in the biochemical reactions and the 
parameter uncertainty of the models are important ques- 
tions that would be explored. Furthermore, the design of 
circuits could be combined with tools for the design of 
synthetic DNA sequences. This would exploit the inter- 
actions between nucleic acids and the reengineering of 
natural proteins. Promoters with targeted transcription 
rates or multiple operators (35,36), small RNAs with 
targeted secondary structures (6), or chimeric proteins 
acting as new transcription factors (74) are examples of 
what we could design computationally. These all elements 
would be modeled by transfer functions and these would 
be stored in a library (virtual or real). Importantly, it 
could be also specified a degree of evolvability, by which 
the value of the kinetic parameters characterizing that 
element would be susceptible to be changed after specific 
mutations in its sequence. Finally, the cellular chassis in 
which the circuit is going to be deployed could be also 
introduced as a generalized element by modeling the 
host elements that require the circuit for its expression 
(78). This would allow to provide a prediction of the 
response of the engineered cell under the conditions for 
which the circuit was designed, and consequently improve 
the design process. 
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